Has slots for bbox which are the spatial extents, and a
CRS object for the coordinate reference system.
getClass("Spatial")
## Class "Spatial" [package "sp"]
##
## Slots:
##
## Name: bbox proj4string
## Class: matrix CRS
##
## Known Subclasses:
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2
getClass("CRS")
## Class "CRS" [package "sp"]
##
## Slots:
##
## Name: projargs
## Class: character
# Example instantiation
bb <- matrix(c(30, 85, 300, 90), ncol = 2, dimnames = list(NULL,c("min", "max")))
crs <- CRS("+proj=longlat")
Spatial(bb, crs)
## class : Spatial
## features : 1
## extent : 30, 300, 85, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
Class SpatialPoints adds a coords slot to
Spatial
getClass("SpatialPoints")
## Class "SpatialPoints" [package "sp"]
##
## Slots:
##
## Name: coords bbox proj4string
## Class: matrix matrix CRS
##
## Extends: "Spatial", "SpatialVector"
##
## Known Subclasses:
## Class "SpatialPointsDataFrame", directly
## Class "SpatialPixels", directly
## Class "SpatialPixelsDataFrame", by class "SpatialPixels", distance 2
# example instantiation
cran <- read.table("data/CRAN051001a.txt", header = TRUE)
cran_mat <- cbind(cran$long, cran$lat)
llcrs <- CRS("+proj=longlat +ellps=WGS84")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
## definition
cran_sp <- SpatialPoints(cran_mat, proj4string = llcrs)
summary(cran_sp)
## Object of class SpatialPoints
## Coordinates:
## min max
## coords.x1 -122.95000 153.0333
## coords.x2 -37.81667 57.0500
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
# accessor methods
coordinates(cran_sp)[which(cran$loc == "Brazil"),] # index coordinates
## coords.x1 coords.x2
## [1,] -49.26667 -25.41667
## [2,] -42.86667 -20.75000
## [3,] -43.20000 -22.90000
## [4,] -47.63333 -22.71667
## [5,] -46.63333 -23.53333
bbox(cran_sp) # bounding box
## min max
## coords.x1 -122.95000 153.0333
## coords.x2 -37.81667 57.0500
summary(cran_sp[which(cran$loc == "Brazil"), ]) # can index class directly, new bbox.
## Object of class SpatialPoints
## Coordinates:
## min max
## coords.x1 -49.26667 -42.86667
## coords.x2 -25.41667 -20.75000
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 5
Used for matching coordinates with other covariates, like name of coordinate. Needs to have matching row names. Objects are meant to operate like dataframes.
# matches by rownames
row.names(cran_mat) <- 1:nrow(cran)
# cran %>% row.names() # matches
cran_spdf <- SpatialPointsDataFrame(cran_mat, cran, proj4string = llcrs, match.ID = TRUE)
summary(cran_spdf) # meant to
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -122.95000 153.0333
## coords.x2 -37.81667 57.0500
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
## Data attributes:
## place north east loc
## Length:54 Length:54 Length:54 Length:54
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## long lat
## Min. :-122.9500 Min. :-37.82
## 1st Qu.: -47.3833 1st Qu.: 34.52
## Median : 7.8500 Median : 42.73
## Mean : -0.6617 Mean : 31.73
## 3rd Qu.: 16.8333 3rd Qu.: 47.58
## Max. : 153.0333 Max. : 57.05
# alternate way, just attach coordinates and proj4string by assignment
cran_spdf_alt <- cran
# coordinates(cran_spdf_alt) <- c("long", "lat") # can select out columns if already in dataframe
coordinates(cran_spdf_alt) <- cran_mat # coordinates as matrix
proj4string(cran_spdf_alt) <- llcrs # crs
summary(cran_spdf_alt) # Turns it into a SpatialPointsDataFrame by assigning coordinate/proj
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -122.95000 153.0333
## coords.x2 -37.81667 57.0500
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
## Data attributes:
## place north east loc
## Length:54 Length:54 Length:54 Length:54
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## long lat
## Min. :-122.9500 Min. :-37.82
## 1st Qu.: -47.3833 1st Qu.: 34.52
## Median : 7.8500 Median : 42.73
## Mean : -0.6617 Mean : 31.73
## 3rd Qu.: 16.8333 3rd Qu.: 47.58
## Max. : 153.0333 Max. : 57.05
The spacetime package provides dataframes for certain types of spatial data.
STFDFgetClass("ST")
## Class "ST" [package "spacetime"]
##
## Slots:
##
## Name: sp time endTime
## Class: Spatial xts POSIXct
##
## Known Subclasses:
## Class "STF", directly
## Class "STS", directly
## Class "STI", directly
## Class "STT", directly
## Class "STFDF", by class "STF", distance 2
## Class "STSDF", by class "STS", distance 2
## Class "STIDF", by class "STI", distance 2
## Class "STTDF", by class "STT", distance 2
# Inheritance
This section follows the tutorial from Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny
data("scotland")
# getClass("SpatialPolygons") # has polygons, plotOrder, bbox, proj4string
# slotNames(map) # getting the value names of the object itself
map <- scotland$spatial.polygon
plot(map) # class SpatialPolygons
The map is in the projection OSGB 1936/British National Grid which has EPSG code 27700. The proj4 string of this projection can be seen in https://spatialreference.org/ref/epsg/27700/proj4/ or can be obtained with R as follows:
codes <- rgdal::make_EPSG()
codes[which(codes$code == "27700"), ]
## code note
## 6288 27700 OSGB36 / British National Grid
## prj4
## 6288 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
## prj_method
## 6288 Transverse Mercator
# assign
proj4string(map) <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=km +no_defs"
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
## Warning in ReplProj4string(obj, CRS(value)): A new CRS was assigned to an object with an existing CRS:
## +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform
# reproject into WGS84 for leaflet
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84 +no_defs"))
Data preparation
county: id of each countryY: observed number of lip cancer cases in each
countyE: expected number of lip cancer cases in each
countyAFF: proportion of population in agriculture, fishing
and forestrySIR: standardized incidence ratio# prepare the dataset and calculate SIR
d <- scotland$data[,c("county.names", "cases", "expected", "AFF")]
names(d) <- c("county", "Y", "E", "AFF")
d$SIR <- d$Y / d$E # observed over expected
# prepare SpatialPolygonsDataFrame, SpatialPolygons + Data basically
rownames(d) <- d$county # set rownames to county for matching later
# getClass("SpatialPolygonsDataFrame") # just adds "data" slot to SpatialPolygons
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE) # make Object for leaflet
l <- leaflet(map) %>% addTiles()
# domain gives possible values to be mapped, returns function
pal <- colorNumeric(palette = "YlOrRd", domain = map$SIR)
# Create the tooltip label
labels <- sprintf("<strong> %s </strong> <br/>
Observed: %s <br/> Expected: %s <br/>
AFF: %s <br/> SIR: %s",
map$county, map$Y, round(map$E, 2),
map$AFF, round(map$SIR, 2)
) %>%
lapply(htmltools::HTML)
l %>%
addPolygons(
color = "grey", weight = 1, # the outlines
fillColor = ~ pal(SIR), fillOpacity = 0.5, # how to fill the shapes
highlightOptions = highlightOptions(weight = 4), # thicker outline on hover
label = labels, # the tooltip label
labelOptions = labelOptions( # some general options for the labels
style = list(
"font-weight" = "normal",
padding = "3px 8px"
),
textsize = "15px", direction = "auto"
)
) %>%
# legend in bottom right
addLegend(
pal = pal, values = ~SIR, opacity = 0.5,
title = "SIR", position = "bottomright"
)
For the modeling piece, we model with CAR,
\[ \begin{aligned} Y_i &\sim Poisson(E_i \theta_i), i = 1,\dots , n \\ \log(\theta_i) &= \beta_0 + \beta_1 \times AFF_i + u_i + v_i \end{aligned} \] - \(E_i\) is expected count - \(\theta_i\) is relative risk in area \(i\) - \(AFF_i\) is the fishing/farming covariate - \(u_i\) is spatial component with conditional guassian model, we assum \(u_i|\mathbf{u_{-i}} \sim N(\bar u_{\delta_i}, \frac{\sigma^2_u}{n_{\delta_i}})\) - \(v_i\) is unstructured spatial effect
# create neighborhood matrix from scottish data
nb <- poly2nb(map) # "nb" class
head(nb) # just a list of neighbors by node
## [[1]]
## [1] 5 9 19
##
## [[2]]
## [1] 7 10
##
## [[3]]
## [1] 12
##
## [[4]]
## [1] 18 20 28
##
## [[5]]
## [1] 1 12 19
##
## [[6]]
## [1] 0
nb # printing shows some summary statistics of the dataset, as well as "disconnected regions" with no neighbors
## Neighbour list object:
## Number of regions: 56
## Number of nonzero links: 234
## Percentage nonzero weights: 7.461735
## Average number of links: 4.178571
## 3 regions with no links:
## orkney shetland western.isles
# nb2mat(nb, style="B",zero.policy=T) # convert to an adjacency matrix
Can calculate the Moran’s I, allows testing of spatial autocorrelation in the data, defined as
\[ \begin{equation} I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}} \frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i - \bar{x})^2} \end{equation} \]
set.ZeroPolicyOption(TRUE) # setting globally, in argument it didn't work
## [1] FALSE
nb_W <- nb2listw(nb, style = "W")
nb_global <- Szero(nb_W) # get global sum of weights
# direct calculation of the statistic
moran(x =map$SIR,
listw = nb_W,
n = length(map$SIR),
S0 = nb_global)
## $I
## [1] 0.5251582
##
## $K
## [1] 5.023736
moran.plot(map$SIR, nb_W) # upper right and lower left are positively correlated, slope is ~ approximately the moran statistic
moran.mc(map$SIR, nb_W, 100) # Lots of spatial correlation, randomization test.
##
## Monte-Carlo simulation of Moran I
##
## data: map$SIR
## weights: nb_W
## number of simulations + 1: 101
##
## statistic = 0.49702, observed rank = 101, p-value = 0.009901
## alternative hypothesis: greater
Not quite sure how to interpret this plot,
# transfer the neighbors object to INLA
nb2INLA("output/map.adj", nb) # export to file, in INLA format
g <- inla.read.graph(filename = "output/map.adj") # read the .adj file for INLA
# indices for random effects u, v
map$idareau <- 1:nrow(map@data) # u
map$idareav <- 1:nrow(map@data) # v
# f() denotes random effect, "besag"
formula <- Y ~ AFF +
f(idareau, model = "besag", graph = g, scale.model = TRUE) + # spatial dependence
f(idareav, model = "iid") # unstructured noise
# fit the formual
res <- inla(formula, # formula from above
family = "poisson", # log link
data = map@data, # the data file
E = E, # "Known component in the mean for the Poisson likelihoods defined as E exp(eta)" see docs
control.predictor = list(compute = TRUE), # compute posteriors of predictions
control.compute=list(return.marginals.predictor=TRUE) # marginals of each predictor as well, for exceedence probabilities
)
summary(res) # output
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 3.52, Running = 5.56, Post = 0.0319, Total = 9.12
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.305 0.120 -0.538 -0.305 -0.068 NA 0
## AFF 4.327 1.276 1.746 4.353 6.765 NA 0
##
## Random effects:
## Name Model
## idareau Besags ICAR model
## idareav IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for idareau 4.14 1.44 2.02 3.91 7.62 NA
## Precision for idareav 19639.46 19918.35 1375.78 13644.97 72579.48 NA
##
## Marginal log-Likelihood: -189.72
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
The results are pretty standard for interpretation, we can get
marginals of the parameters (AFF) with smoothed (spline)
marginals
# results has an output of marginals for each of the fixed effects, on some grid
# plot(res$marginals.fixed$AFF)
# lines(marginal$x, marginal$y)
marginal <- inla.smarginal(res$marginals.fixed$AFF) # smoothed marginals
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() +
labs(x = expression(beta[1]), y = "Density") +
geom_vline(xintercept = 0, col = "black") + theme_bw()
Each fitted value for each of the counties can be extracted from the summaries.
head(res$summary.fitted.values)
## mean sd 0.025quant 0.5quant 0.975quant mode
## fitted.Predictor.01 4.965762 1.4504631 2.643878 4.790508 8.291851 NA
## fitted.Predictor.02 4.396906 0.6752405 3.172452 4.362819 5.817591 NA
## fitted.Predictor.03 3.620712 1.0169095 1.918926 3.523866 5.883690 NA
## fitted.Predictor.04 3.083339 0.8946611 1.627708 2.983424 5.113472 NA
## fitted.Predictor.05 3.329251 0.7500895 2.042165 3.267024 4.974026 NA
## fitted.Predictor.06 2.975147 0.9194491 1.490216 2.869251 5.070599 NA
final version of the leaflet map, with relative risk, and 95% CI for RR.
# Add to map df data
map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]
# Leaflet code
pal <- colorNumeric(palette = "YlOrRd", domain = map$RR)
labels <- sprintf("<strong> %s </strong> <br/>
Observed: %s <br/> Expected: %s <br/>
AFF: %s <br/> SIR: %s <br/> RR: %s (%s, %s)",
map$county, map$Y, round(map$E, 2),
map$AFF, round(map$SIR, 2), round(map$RR, 2),
round(map$LL, 2), round(map$UL, 2)
) %>% lapply(htmltools::HTML)
lRR <- leaflet(map) %>%
addTiles() %>%
addPolygons(
color = "grey", weight = 1, fillColor = ~ pal(RR),
fillOpacity = 0.5,
highlightOptions = highlightOptions(weight = 4),
label = labels,
labelOptions = labelOptions(
style =
list(
"font-weight" = "normal",
padding = "3px 8px"
),
textsize = "15px", direction = "auto"
)
) %>%
addLegend(
pal = pal, values = ~RR, opacity = 0.5, title = "RR",
position = "bottomright"
)
lRR
# exceedence probabilities, cut offs for action
# P(\theta_i > c), probability that a specific area i has relative risk higher than some threshold.
sapply(res$marginals.fitted.values, FUN = function(marg){1 - inla.pmarginal(q = 2, marginal = marg)})
## fitted.Predictor.01 fitted.Predictor.02 fitted.Predictor.03 fitted.Predictor.04
## 9.975709e-01 9.999996e-01 9.667358e-01 9.067867e-01
## fitted.Predictor.05 fitted.Predictor.06 fitted.Predictor.07 fitted.Predictor.08
## 9.793448e-01 8.680805e-01 9.780170e-01 4.596383e-01
## fitted.Predictor.09 fitted.Predictor.10 fitted.Predictor.11 fitted.Predictor.12
## 5.130059e-01 9.720600e-01 6.550928e-01 9.442041e-01
## fitted.Predictor.13 fitted.Predictor.14 fitted.Predictor.15 fitted.Predictor.16
## 7.405638e-01 5.351114e-01 4.747477e-01 6.258552e-01
## fitted.Predictor.17 fitted.Predictor.18 fitted.Predictor.19 fitted.Predictor.20
## 5.192967e-01 1.280056e-02 4.852896e-01 8.802111e-02
## fitted.Predictor.21 fitted.Predictor.22 fitted.Predictor.23 fitted.Predictor.24
## 3.173132e-02 2.622425e-02 1.326642e-02 6.028813e-05
## fitted.Predictor.25 fitted.Predictor.26 fitted.Predictor.27 fitted.Predictor.28
## 3.415062e-03 1.684688e-03 1.839832e-03 2.110785e-03
## fitted.Predictor.29 fitted.Predictor.30 fitted.Predictor.31 fitted.Predictor.32
## 1.553052e-03 5.385300e-05 2.241450e-04 2.026798e-01
## fitted.Predictor.33 fitted.Predictor.34 fitted.Predictor.35 fitted.Predictor.36
## 1.940328e-03 0.000000e+00 3.589546e-05 1.326963e-04
## fitted.Predictor.37 fitted.Predictor.38 fitted.Predictor.39 fitted.Predictor.40
## 4.311477e-05 0.000000e+00 5.081891e-03 4.865651e-06
## fitted.Predictor.41 fitted.Predictor.42 fitted.Predictor.43 fitted.Predictor.44
## 0.000000e+00 0.000000e+00 6.511447e-03 0.000000e+00
## fitted.Predictor.45 fitted.Predictor.46 fitted.Predictor.47 fitted.Predictor.48
## 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## fitted.Predictor.49 fitted.Predictor.50 fitted.Predictor.51 fitted.Predictor.52
## 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## fitted.Predictor.53 fitted.Predictor.54 fitted.Predictor.55 fitted.Predictor.56
## 0.000000e+00 0.000000e+00 2.773613e-03 5.010538e-04
Had to pull the data directly from the archived package, it seems to have been taken off CRAN because the contact author is no longer active. Link to Archive
dohio <- read.csv("data/dataohiocomplete.csv")
head(dohio)
## county gender race year y n NAME
## 1 1 1 1 1968 6 8912 Adams
## 2 1 1 1 1969 5 9139 Adams
## 3 1 1 1 1970 8 9455 Adams
## 4 1 1 1 1971 5 9876 Adams
## 5 1 1 1 1972 8 10281 Adams
## 6 1 1 1 1973 5 10876 Adams
map <- as_Spatial(st_read("data/fe_2007_39_county", layer="fe_2007_39_county")) # read shape files, not quite, then convert to Spatial for sp methods
## Reading layer `fe_2007_39_county' from data source
## `/Users/mliou/Box Sync/gradwisc/reference/cheatsheets/mydocs/website/static/rmd/spatial/data/fe_2007_39_county'
## using driver `ESRI Shapefile'
## Simple feature collection with 88 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -84.8203 ymin: 38.40342 xmax: -80.5182 ymax: 42.32713
## Geodetic CRS: NAD83
plot(map)
# observed cases
d <- aggregate(
x = dohio$y,
by = list(county = dohio$NAME, year = dohio$year),
FUN = sum
)
names(d) <- c("county", "year", "Y")
head(d)
## county year Y
## 1 Adams 1968 6
## 2 Allen 1968 32
## 3 Ashland 1968 15
## 4 Ashtabula 1968 27
## 5 Athens 1968 12
## 6 Auglaize 1968 7
# reorder correctly,
dohio <- dohio[order(
dohio$county,
dohio$year,
dohio$gender,
dohio$race
), ]
# check sorting
# dohio[1:20, ]
# 2 race x 2 gender
n.strata <- 4
E <- expected( # function needs sorted data
population = dohio$n,
cases = dohio$y,
n.strata = n.strata
)
nyears <- length(unique(dohio$year)) # number of years
countiesE <- rep(unique(dohio$NAME), # county vector
each = nyears)
ncounties <- length(unique(dohio$NAME))
yearsE <- rep(unique(dohio$year), # rep by counties
times = ncounties)
# putting it all together
dE <- data.frame(county = countiesE,
year = yearsE,
E = E)
head(dE)
## county year E
## 1 Adams 1968 8.278660
## 2 Adams 1969 8.501767
## 3 Adams 1970 8.779313
## 4 Adams 1971 9.175276
## 5 Adams 1972 9.548736
## 6 Adams 1973 10.099777
# put it back into the main dataset
d <- merge(d, dE, by = c("county", "year"))
head(d)
## county year Y E
## 1 Adams 1968 6 8.278660
## 2 Adams 1969 5 8.501767
## 3 Adams 1970 9 8.779313
## 4 Adams 1971 6 9.175276
## 5 Adams 1972 10 9.548736
## 6 Adams 1973 7 10.099777
# calc SIR
d$SIR <- d$Y / d$E
head(d)
## county year Y E SIR
## 1 Adams 1968 6 8.278660 0.7247549
## 2 Adams 1969 5 8.501767 0.5881130
## 3 Adams 1970 9 8.779313 1.0251372
## 4 Adams 1971 6 9.175276 0.6539313
## 5 Adams 1972 10 9.548736 1.0472591
## 6 Adams 1973 7 10.099777 0.6930846
# reshape into lengthwise file with rows as counties and columns as years
dw <- reshape(d,
timevar = "year",
idvar = "county",
direction = "wide"
)
map <- merge(map, dw, by.x = "NAME", by.y = "county")
map_sf <- st_as_sf(map) # into sf for ggplot
map_sf <- gather(map_sf, year, SIR, paste0("SIR.", 1968:1988)) # longwise data now
map_sf$year <- as.integer(substring(map_sf$year, 5, 8)) # just get the year string
# make plot
ggplot(map_sf) +
geom_sf(aes(fill = SIR)) + # map the shapefiles
facet_wrap(~year, dir = "h", ncol = 7) +
ggtitle("SIR") + theme_bw() +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank()
) +
scale_fill_gradient2(
midpoint = 1, low = "blue", mid = "white", high = "red"
)
ggplot(d, aes(x = year, y = SIR,
group = county, color = county)) +
geom_line() + geom_point(size = 2) + theme_bw() +
theme(legend.position = "none") +
gghighlight(county == "Adams")
## Warning: Tried to calculate with group_by(), but the calculation failed.
## Falling back to ungrouped filter operation...
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## Warning: Tried to calculate with group_by(), but the calculation failed.
## Falling back to ungrouped filter operation...
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## label_key: county
Now we actually model the data with a spatio-temporal component
nb <- poly2nb(map) # convert to adjmat
nb2INLA("data/map_ohio.adj", nb)
g <- inla.read.graph(filename = "data/map_ohio.adj")
d$idarea <- as.numeric(as.factor(d$county))
d$idarea1 <- d$idarea # make copy of the idarea, becuase two terms with spatial random componenet, but in INLA, only one variable can be associated to the `f()` at a time.
d$idtime <- 1 + d$year - min(d$year) # index years/dates from 1
# for documentation of how this matches with model
# inla.doc(what = "bym")
formula <- Y ~ f(idarea, model = "bym", graph = g) + # spatial component
f(idarea1, idtime, model = "iid") + # time x area interaction
idtime # global trend across time
res <- inla(formula,
family = "poisson",
data = d,
E = E,
control.predictor = list(compute = TRUE)
)
summary(res)
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 4.1, Running = 1.91, Post = 0.0435, Total = 6.05
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.530 0.024 -0.578 -0.530 -0.483 NA 0
## idtime 0.037 0.001 0.035 0.037 0.039 NA 0
##
## Random effects:
## Name Model
## idarea BYM model
## idarea1 IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for idarea (iid component) 24.53 4.21 17.20 24.21
## Precision for idarea (spatial component) 1890.38 1772.12 121.01 1358.39
## Precision for idarea1 49638.84 17544.85 24285.89 46629.35
## 0.975quant mode
## Precision for idarea (iid component) 33.75 NA
## Precision for idarea (spatial component) 6552.57 NA
## Precision for idarea1 92463.65 NA
##
## Marginal log-Likelihood: -5940.95
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
It’s hard to interpret these
# plot with RR values
d$RR <- res$summary.fitted.values[, "mean"]
d$LL <- res$summary.fitted.values[, "0.025quant"]
d$UL <- res$summary.fitted.values[, "0.975quant"]
map_sf <- merge(
map_sf, d,
by.x = c("NAME", "year"),
by.y = c("county", "year")
)
ggplot(map_sf) + geom_sf(aes(fill = RR)) +
facet_wrap(~year, dir = "h", ncol = 7) +
ggtitle("RR") + theme_bw() +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank()
) +
scale_fill_gradient2(
midpoint = 1, low = "blue", mid = "white", high = "red"
)
This exposition follows Spatio-Temporal Statistics With R (Wikle, Zammit-Mangion, Cressie)
This example assumes the Itero-Differential Equation (IDE) is given, and we will show how to simulate from a given IDE. The ID we consider is
\[ \begin{aligned} Y_t(s) = \int_{D_s} m(s, x; \mathbf{\theta_p})Y_{t-1}(x)\,dx + \eta_t(s), \quad s,x\in D_s \end{aligned} \]
set.seed(1)
# create grid for simulation, (between 0,1)
ds <- 0.01 # step for spatial grid
s_grid <- seq(0, 1, by = ds) # spatial grid
N <- length(s_grid)
nT <- 201 # number of time points
t_grid <- 0:(nT-1) # time grid
st_grid <- expand.grid(s = s_grid, t = t_grid) # space-time grid
m <- function(x, s,thetap) {
gamma <- thetap[1] # amplitude
l <- thetap[2] # length
offset <- thetap[3] # offset
D <- outer(s + offset, x, "-") # displacements
gamma * exp(-D^2/l) # kernel eval.
}
# see how the parameterization differs for a few
thetap <- list()
thetap[[1]] <- c(40, 0.0002, 0) # narrow, centered.
thetap[[2]] <- c(5.75, 0.01, 0) # wider
thetap[[3]] <- c(8, 0.005, 0.1) # shifted right
thetap[[4]] <- c(8, 0.005, -0.1) # shifted left
m_x <- m(s=0.5, x = s_grid,
thetap=thetap[[1]]) %>%
as.numeric()
df <- data.frame(x=s_grid, m = m_x)
# plot the data
ggplot(df) +
geom_line(aes(x,m)) +
theme_bw()
# The covariance of the spatial noise
# Simulate from Multivariate Normal, by cholesky decomp of the covariance
Sigma_eta <- 0.1 * exp( -abs(outer(s_grid, s_grid, '-') / 0.1))
L <- t(chol(Sigma_eta)) # chol() returns upper Cholesky factor
sim <- L %*% rnorm(nrow(Sigma_eta)) # simulate realization of eta
qplot(s_grid, sim, geom = "line") +
theme_minimal() +
labs(title = "realization of a random spatial process in 1d")
Y <- list()
# outer loop does 4 kernels
for(i in 1:4) {
M <- m(s_grid, s_grid, thetap[[i]]) # create kernel function
Y[[i]] <- data.frame(s = s_grid, # initialize grid and process
t = 0,
Y = 0)
# for each timepoint
for (j in t_grid[-1]) {
prev_Y <- filter(Y[[i]], t == j - 1)$Y # get Y at t-1
eta <- L %*% rnorm(N) # simulate eta
new_Y <- (M %*% prev_Y * ds + eta) %>% as.numeric() # euler approximation, approximating integral with sum
Y[[i]] <- rbind(Y[[i]],
data.frame(s = s_grid,
t = j,
Y = new_Y))
}
}
hovmoller_plot <- function(idx) {
ggplot(Y[[idx]]) +
geom_tile(aes(s,t,fill = Y)) +
scale_x_continuous(expand = c(0,0)) +
scale_y_reverse(expand = c(0,0)) +
fill_scale(name = "Y") +
theme_bw() +
theme(legend.position = "bottom")
}
p1 <- hovmoller_plot(1)
p2 <- hovmoller_plot(2)
p3 <- hovmoller_plot(3)
p4 <- hovmoller_plot(4)
plot_grid(p1, p2, p3, p4, nrow = 1)
Now we learn how to simulate data. The above graphs show the entire realization of the process model. In reality, we don’t get to observe all that data, instead, we just get little bits of the dynamics. So we create an “Incidence Matrix” \(\mathbf{H}_t\) to account for the data that we actually do observe. i.e., the above was the “process model”, now we implement the “data model”. We also assume we make some measurement errors, that are iid w/ mean zero.
\[ \begin{aligned} Z_t = H_tY_t + \varepsilon_t \end{aligned} \]
nobs <- 50 # number of observations at each timepoint
sobs <- sample(s_grid, nobs) # sample points on the grid
# should be w/ replacment right?
# Initialize incidence, nobs x N matrix
Ht <- matrix(0, nrow = nobs, ncol = N) # the incidence matrix
for(i in 1:nobs){
idx <- which(sobs[i] == s_grid) # find elements and set to 1
Ht[i, idx] <- 1
}
z_df <- data.frame() # init data frame to save noisy observations
# for each time point, we get the underlying process, take our observations and add measurement error
for (j in 0:(nT-1)) {
Yt <- Y[[1]] %>% filter(t == j) %>% pull(Y) # get simulated process
zt <- Ht %*% Yt + rnorm(nobs) # map
z_df <- rbind(z_df, data.frame(s = sobs, t = j, z = zt))
}
ggplot(z_df) + geom_point(aes(s, t, colour = z)) +
col_scale(name = "z") +
scale_y_reverse(expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
theme_test() +
labs(title = "Noisy Observation of Latent Process")
First example if simulating from a spatially invariant kernel, this
is what k_spat_invariant=1 specificies. The movement here
is very simple, it simply drifts away.
SIM1 <- simIDE(T = 10, nobs = 100, k_spat_invariant = 1)
str(SIM1, max.level=1)
## List of 6
## $ s_df :'data.frame': 16810 obs. of 4 variables:
## $ z_df :'data.frame': 1000 obs. of 4 variables:
## $ z_STIDF :Formal class 'STIDF' [package "spacetime"] with 4 slots
## $ g_truth :List of 9
## ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
## $ g_obs :List of 9
## ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
## $ IDEmodel:List of 5
## ..- attr(*, "class")= chr "IDE"
s_df - is the space dfz_df - is the observed noisy datag_truth - a ggplot of the true processg_obs - a ggplot of the observationsprint(SIM1$g_truth)
print(SIM1$g_obs)
# returns object of class IDE and initial parameter estimates of \alpha_t
IDEmodel <- IDE(f=z~s1+s2,
data = SIM1$z_STIDF,
dt = as.difftime(1, units = "days"),
grid_size = 41)
# Fits with deOptim, a derivative free optimzation method (same category as Nelder Mead and Simulated Annealing, but with different principles)
# Stands for differential evolution
# function is incredibly intenstive for calculation so the results can also be loaded below -
# fit_results_sim1 <- fit.IDE(IDEmodel, parallelType = 1)
data("IDE_Sim1_results", package = "STRbook")
show_kernel(fit_results_sim1$IDEmodel) # show the kernel of the function
## Kernel is spatially invariant, plotting it centred on the origin.
# parameters of the model can be extracted from the object here
# real values are c(150, .002, -0.1, 0.1)
fit_results_sim1$IDEmodel$get("k") %>% unlist()
## par1 par2 par3 par4
## 152.836345912 0.001977115 -0.101601099 0.100368743
coef(fit_results_sim1$IDEmodel) # true values are c(.2, .2, .2)
## Intercept s1 s2
## 0.2073442 0.1966224 0.1907062
Now we extract the eigenvalues of the fitted evolution matrix to understand its long term behavior.
abs_ev <- eigen(fit_results_sim1$IDEmodel$get("M"))$values %>%
abs()
abs_ev %>% head()
## [1] 0.4640858 0.4640858 0.4492270 0.4492270 0.4289470 0.4289470
The largest of these values is greater than 1, thus we conclude that the dynamics of this evolution are stable.
ST_grid_df <- predict(fit_results_sim1$IDEmodel) # can also take a new prediction grid, but by default it takes the simulation grid that we fit the data model for.
# Prediction map
gpred <- ggplot(ST_grid_df) + # Plot the predictions geom_tile(aes(s1, s2, fill=Ypred)) +
geom_tile(aes(s1,s2, fill = Ypred)) +
facet_wrap(~t) +
fill_scale(name = "Ypred", limits = c(-0.1, 1.4)) +
coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
labs(title = "Predictions")
# Prediction standard error map
gpredse <- ggplot(ST_grid_df) + # Plot the prediction s.es
geom_tile(aes(s1, s2, fill = Ypredse)) +
facet_wrap(~t) +
fill_scale(name = "Ypredse") +
coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
labs(title = "Predictions SE")
plot_grid(gpred, gpredse)
### Now we run an example of spatially varying kernel
spatially varying kerneles are necessary in order to model forms of
advection. specified with k_spat_invariant = 0
SIM2 <- simIDE(T = 15, nobs = 1000, k_spat_invariant = 0)
# the two ggplots of the truth and the simulated observations
print(SIM2$g_truth)
print(SIM2$g_obs)
appears to counter clockwise rotate, and then come to standstill in
bottom of the domain. The varying advection that generated this field is
visualized as follows…
show_kernel(SIM2$IDEmodel, scale = 0.2) # scale says change the arrow sizes.
## Kernel is spatially variant, plotting displacements
# bisquare basis functions, same class as used by FRK.
mbasis_1 <- auto_basis(manifold = plane(), #fns, on the plane
data = SIM2$z_STIDF, # data
nres = 1, # 1 resolution
type = 'bisquare')
show_basis(mbasis_1)
## Note: show_basis assumes spherical distance functions when plotting
# now in our kernel, we can assume that amplitude (thetam1) and the scale (thetam2) are both constant
kernel_basis <- list(thetam1 = constant_basis(),
thetam2 = constant_basis(),
thetam3 = mbasis_1,
thetam4 = mbasis_1)
# specify the model
IDEmodel <- IDE(f = z ~ s1 + s2 + 1,
data = SIM2$z_STIDF,
dt = as.difftime(1, units = "days"),
grid_size = 41,
kernel_basis = kernel_basis)
# fit the IDE, computationally intensive so load from cache
# fit_results_sim2 <- fit.IDE(IDEmodel, parallelType = 1,
# itermax = 400)
data("IDE_Sim2_results", package = "STRbook")
show_kernel(fit_results_sim2$IDEmodel)
## Kernel is spatially variant, plotting displacements
ST_grid_df2 <- predict(fit_results_sim2$IDEmodel)
gpred <- ggplot(ST_grid_df2) + # Plot the predictions geom_tile(aes(s1, s2, fill=Ypred)) +
geom_tile(aes(s1, s2, fill = Ypred)) +
facet_wrap(~t) +
fill_scale(name = "Ypred", limits = c(-0.1, 1.4)) +
coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
labs(title = "Predictions")
These are the predictions that were
Above we did a simulation example, but here we do a real data example
data("radar_STIDF", package = "STRbook")
hindcast how many time intervals in to past to
estimateforecase how many time intervals into the future to
estimate.IDEmodel <- IDE(f = z~1,
data = radar_STIDF,
dt=as.difftime(10, units = "mins"),
grid_size = 41,
forecast = 2,
hindcast = 2)
# fitmodel
# fit_results_radar <- fit.IDE(IDEmodel, parallelType = 1)
data("IDE_Radar_results", package = "STRbook") # load cached results
show_kernel(fit_results_radar$IDEmodel)
## Kernel is spatially invariant, plotting it centred on the origin.
shifted off center, so suggestive of predominatly north east direction.
shift_pars <- (fit_results_radar$IDEmodel$get("k") %>% unlist())[3:4]
print(shift_pars)
## par3 par4
## -5.488385 -1.860784
We can make a calculation here, \(\sqrt{5.5^2 + 1.9^2} = 5.82\) in 10 minutes, 34.91 km per hour.
abs_ev <- eigen(fit_results_radar$IDEmodel$get("M"))$values %>% abs()
summary(abs_ev)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.008377 0.581449 0.671552 0.615239 0.719497 0.786783
maximum eigenvalue less than 1 means that its stable, but more persistent than the simulation example.
ST_grid_df <- predict(fit_results_radar$IDEmodel)
radar_df$time <- format(radar_df$t, "%H:%M")
ST_grid_df$time <- format(ST_grid_df$t, "%H:%M")
## Add time records with missing data
radar_df <- rbind.fill(radar_df, data.frame(time = c("08:05", "08:15",
"10:25", "10:35")))
## Plot of data, with color scale capped to (-20, 60)
gobs <- ggplot(radar_df) +
geom_tile(aes(s1, s2, fill = pmin(pmax(z, -20), 60))) +
fill_scale(limits = c(-20, 60), name = "Z") +
facet_wrap(~time) +
coord_fixed() +
theme_bw() +
labs(title = "Observed Sydney Data")
## Plot of predictions with color scale forced to (-20, 60)
gpred <- ggplot(ST_grid_df) +
geom_tile(aes(s1, s2, fill = Ypred)) + facet_wrap(~time) + coord_fixed() + theme_bw() + fill_scale(limits = c(-20, 60), name = "Ypred") +
labs(title = "Predicted Data, with 2 time periods before and after")
plot_grid(gobs, gpred)
## Warning: Removed 4 rows containing missing values (geom_tile).
We use EOFs are ideal basis functions for exploring the observed signal.
We fit two models, a traditional one in time-series framework with vector autoregression and use method of moments.
# trucate at April 1997 in order to forcast 6 months ahead
data("SSTlandmask", package = "STRbook")
data("SSTlonlat", package = "STRbook")
data("SSTdata", package = "STRbook")
delete_rows <- which(SSTlandmask == 1) # remove land values
SST_Oct97 <- SSTdata[-delete_rows, 334] # save Oct 1997 SSTs
SSTdata <- SSTdata[-delete_rows, 1:328] # until April 1997 , 399 times points normally, but truncate
SSTlonlat$mask <- SSTlandmask # assign mask to df
# 2520 x 328 dataframe
SSTdata %>% head()
## V1 V2 V3 V4 V5 V6 V7
## 16 0.2296104 0.2515602 0.4510937 0.2989082 -0.28913879 -0.20289230 -0.1392193
## 17 0.4862499 0.4507828 0.5332813 0.5499210 -0.22593880 -0.05968857 0.2024994
## 18 0.7855454 0.6135120 0.5860958 0.6160946 -0.11890602 -0.01961136 0.3596077
## 19 0.9411717 0.7603111 0.7017193 0.5956268 -0.02531242 -0.02273369 0.4096889
## 20 1.0239792 0.7840633 0.7537537 0.6158600 0.02359390 -0.06570435 0.3028908
## 21 1.0313301 0.8165626 0.8502331 0.6375790 0.09546661 -0.03921890 0.2607803
## V8 V9 V10 V11 V12 V13
## 16 -0.3317967 -0.08953095 0.000623703 0.16234589 0.3113289 0.6296101
## 17 0.0696888 0.15523529 0.237970350 0.11281204 0.5346890 0.7362499
## 18 0.2849217 0.20914078 0.395938870 -0.06546974 0.5482826 0.7855454
## 19 0.3793736 0.20929337 0.427658080 -0.05945587 0.5828152 0.8411713
## 20 0.4363270 0.26359367 0.372110370 -0.05234337 0.5239067 0.9239807
## 21 0.4079704 0.33273506 0.348985670 -0.01726532 0.4121094 0.9813290
## V14 V15 V16 V17 V18 V19
## 16 -0.1484413 -0.29890633 -0.4010925 -0.18914032 0.1471081 -0.1392174
## 17 0.2507839 -0.16671944 -0.1000786 0.02406311 0.5903110 0.3024979
## 18 0.7135143 -0.06390381 0.1160946 0.18109322 0.9303894 0.5096092
## 19 1.0103111 0.05171967 0.2456264 0.32468605 1.0772667 0.5596886
## 20 1.0840626 0.10375214 0.3158608 0.47359467 1.0842972 0.5528908
## 21 1.0665646 0.20023537 0.3875771 0.69546700 1.1107807 0.5607834
## V20 V21 V22 V23 V24 V25
## 16 -0.08179855 -0.2395306 -0.1993771 -0.58765602 -0.4886703 -0.3203907
## 17 0.06968880 0.1552353 0.1879711 -0.23718643 -0.3653126 -0.2137508
## 18 0.18492317 0.4591408 0.4959373 -0.01547051 -0.3517170 -0.3144531
## 19 0.22937393 0.6592922 0.6776581 0.14054489 -0.3171864 -0.5088291
## 20 0.28632736 0.7135944 0.6721096 0.14765739 -0.2760944 -0.6260185
## 21 0.30796814 0.6827354 0.6489849 0.08273506 -0.3378906 -0.6186695
## V26 V27 V28 V29 V30 V31
## 16 -0.4984398 -0.5489063 -0.10109329 -0.53913879 -0.3528919 -0.08921814
## 17 -0.3492184 -0.6167183 0.14992332 -0.22593880 -0.2096882 -0.14750099
## 18 -0.3864861 -0.6639042 0.21609306 0.08109283 -0.1196098 -0.24039078
## 19 -0.4396877 -0.4982815 0.14562798 0.12468910 -0.1227322 -0.24031067
## 20 -0.6159363 -0.2462463 0.06586075 0.12359619 -0.2157040 -0.09710693
## 21 -0.6834374 0.1502342 0.08757782 0.09546661 -0.2392197 0.11078262
## V32 V33 V34 V35 V36 V37 V38
## 16 0.06820297 -0.33953094 0.2006245 0.7123470 -0.03866959 0.3796101 0.7515602
## 17 0.21968651 -0.14476395 0.5379696 0.7128124 -0.31531143 0.5362492 0.9007816
## 18 0.13492393 -0.09085846 0.6959381 0.6345291 -0.50171661 0.6355476 0.9135132
## 19 0.07937241 -0.04070663 0.6276588 0.5905438 -0.41718483 0.6911717 0.8103104
## 20 0.03632736 0.06359291 0.4221096 0.5476570 -0.17609215 0.7239800 0.6840630
## 21 0.05797005 0.33273506 0.2989864 0.6327343 0.06211090 0.6813297 0.5665626
## V39 V40 V41 V42 V43 V44 V45
## 16 0.7010937 -0.4510918 -0.23913956 0.7971077 0.7107811 0.7182026 0.9604683
## 17 0.7832794 -0.2000790 -0.02593803 0.8903122 0.8024979 0.7696876 1.0552330
## 18 0.7860966 0.1160946 0.23109245 0.8303890 0.7096100 0.6849232 0.9091415
## 19 0.7517204 0.3456287 0.52468681 0.6272678 0.6096878 0.6793728 0.7592945
## 20 0.8037529 0.4158592 0.72359467 0.5342960 0.5528908 0.6863270 0.6135922
## 21 0.8002338 0.3875771 0.79546738 0.5607796 0.4607811 0.6579704 0.5327339
## V46 V47 V48 V49 V50 V51 V52
## 16 0.9506245 0.6123447 0.7113304 0.3296108 0.1515617 -0.5489063 -0.001091003
## 17 1.1879711 0.8128128 0.6846886 0.4362488 0.1507816 -0.5167198 0.049921036
## 18 1.2459373 0.9845295 0.5982838 0.4855461 0.1635151 -0.3139057 0.066095352
## 19 1.0776577 0.9905453 0.6328144 0.4911728 0.1603108 -0.1482792 0.145627980
## 20 1.0221100 0.9976559 0.7239056 0.5239792 0.1840630 -0.1462479 0.265859600
## 21 0.9489842 1.0327358 0.8121109 0.5813293 0.2665634 -0.1497650 0.387578960
## V53 V54 V55 V56 V57 V58
## 16 -0.13913918 -0.102891920 -0.089218140 -0.28179550 -0.03952980 -0.3493748
## 17 -0.07593918 -0.009689331 0.002498627 -0.23031044 0.05523491 -0.2120304
## 18 -0.06890678 -0.169610980 0.009609222 -0.21507835 -0.04085922 -0.2540627
## 19 -0.12531281 -0.272733690 0.009689331 -0.12062645 -0.09070587 -0.2223415
## 20 -0.17640495 -0.165702820 0.102891920 -0.11367226 -0.13640785 -0.2778912
## 21 -0.05453491 0.060779572 0.310781480 -0.09202957 -0.21726608 -0.2510128
## V59 V60 V61 V62 V63 V64
## 16 -0.8376541 -1.2386703 -0.77038956 -0.9984398 -0.14890671 -0.20109177
## 17 -0.7871876 -0.9153118 -0.56375122 -0.8492184 0.08328247 0.04992104
## 18 -0.7654705 -0.6017170 -0.41445351 -0.5864868 0.18609428 0.21609497
## 19 -0.6094551 -0.2671871 -0.30882835 -0.4896889 0.25172043 0.09562874
## 20 -0.3523426 0.0239048 -0.22602081 -0.3659363 0.25375366 0.06586075
## 21 -0.1672649 0.3121109 -0.01867104 -0.3334370 0.20023537 0.13757706
## V65 V66 V67 V68 V69 V70 V71
## 16 0.1108608 -0.50289345 0.2607803 0.4182014 0.2604694 -0.34937477 0.212347030
## 17 0.3240623 -0.25968933 0.5024986 0.5196876 0.5052338 -0.21203041 0.212812420
## 18 0.3310947 -0.06961060 0.4596081 0.5849247 0.6091423 -0.05406189 0.134529110
## 19 0.2746868 -0.02273369 0.3096886 0.6793728 0.6092930 0.02765655 0.040544510
## 20 0.2235947 0.03429604 0.3028908 0.6863270 0.5635929 0.07210922 -0.002342224
## 21 0.2454681 0.16077995 0.4107819 0.6079712 0.4327354 0.09898567 0.082735062
## V72 V73 V74 V75 V76 V77
## 16 -0.03866959 -0.17038918 -0.3984394 -0.1989078 -0.1510925 -0.2891388
## 17 0.18468857 -0.06375122 -0.6492176 -0.1167183 -0.3500767 -0.3259392
## 18 0.24828339 -0.01445389 -0.8364849 -0.1139050 -0.6339054 -0.4189053
## 19 0.18281364 -0.05882835 -0.9896908 -0.2482796 -0.8043728 -0.5253124
## 20 0.12390709 -0.12601852 -1.0659370 -0.2962456 -0.7841396 -0.6264038
## 21 0.21210861 0.03132820 -0.9334354 -0.3997669 -0.6624222 -0.6045341
## V78 V79 V80 V81 V82 V83
## 16 0.2971077 -0.23921967 -0.4817963 -0.9395313 0.10062408 -0.1876545
## 17 0.1403122 -0.04750061 -0.4303112 -0.9447651 -0.01202965 -0.1871872
## 18 -0.0696106 0.05960846 -0.5650768 -0.9908581 -0.15406227 -0.2654705
## 19 -0.1727333 0.10968781 -0.6706276 -0.9907055 -0.27234268 -0.4594555
## 20 -0.2157040 0.00289154 -0.5636711 -0.9864063 -0.22789192 -0.6023426
## 21 -0.2392197 -0.33921814 -0.3420296 -0.9172649 -0.15101433 -0.7672653
## V84 V85 V86 V87 V88 V89
## 16 -0.2386704 0.07961082 0.0015602112 -0.1989078 0.19890785 0.01086044
## 17 -0.1653118 -0.06374931 0.0007820129 -0.2167206 0.09992027 -0.07593918
## 18 -0.1017170 -0.26445389 -0.0864849090 -0.3139057 -0.13390541 -0.21890831
## 19 -0.1671848 -0.60882950 -0.2396888700 -0.4982815 -0.30437279 -0.32531166
## 20 -0.2260933 -0.97601891 -0.5159378100 -0.6962471 -0.43413925 -0.47640419
## 21 -0.3878899 -1.11866950 -0.6334362000 -0.8997669 -0.46242142 -0.55453300
## V90 V91 V92 V93 V94 V95
## 16 -0.252893450 -0.6392193 -0.13179779 0.31046867 -0.1993771 -0.6376553
## 17 -0.009689331 -0.5475006 0.06968880 0.15523529 -0.1620293 -0.2371864
## 18 0.180389400 -0.6403923 0.08492279 0.05914116 -0.1540623 0.1345291
## 19 0.227266310 -0.6903114 0.07937241 -0.04070473 -0.3223419 0.3405437
## 20 0.134296420 -0.7971077 0.03632736 -0.13640785 -0.4778900 0.3476582
## 21 0.060779572 -0.9392185 -0.04203033 -0.31726456 -0.5010147 0.1327343
## V96 V97 V98 V99 V100 V101
## 16 -0.4386711 0.12961006 0.35155869 0.001092911 -0.001092911 0.01086044
## 17 -0.4153118 -0.01375008 0.30078316 0.033281326 -0.100078580 -0.07593918
## 18 -0.3517170 -0.21445274 0.11351395 -0.063903809 0.016094208 -0.21890831
## 19 -0.3671856 -0.45882797 0.06031227 -0.298280720 0.045627594 -0.32531166
## 20 -0.4260941 -0.57601929 0.03406334 -0.496246340 -0.034141541 -0.37640381
## 21 -0.4878883 -0.76867104 0.01656151 -0.699766160 -0.112421040 -0.60453415
## V102 V103 V104 V105 V106 V107
## 16 -0.1028919 -0.4892197 -0.3317967 -0.3895321 -0.4493771 -0.4876537
## 17 -0.1096878 -0.2475014 -0.2303123 -0.1447659 -0.6620293 -0.3871880
## 18 -0.2196102 -0.2403908 -0.5150776 -0.2408600 -0.8540630 -0.5154686
## 19 -0.3727322 -0.3903122 -0.7706261 -0.4407063 -1.0223427 -0.6594562
## 20 -0.5157051 -0.7471085 -1.0136719 -0.5364075 -1.1278896 -0.7523422
## 21 -0.5892181 -1.1392174 -1.1420307 -0.5672646 -1.2010136 -0.7672672
## V108 V109 V110 V111 V112 V113
## 16 -0.08866882 -0.07039070 -0.29843903 -0.24890709 -0.2510929 0.06085968
## 17 0.08468819 0.03624916 -0.14921761 -0.11671829 -0.4000778 0.22406006
## 18 0.14828491 -0.01445389 -0.03648567 -0.01390266 -0.4839039 0.23109245
## 19 -0.06718636 -0.10882950 0.01031113 0.05171967 -0.5543728 0.17468834
## 20 -0.37609291 -0.12601852 -0.01593780 0.05375290 -0.4841404 0.17359543
## 21 -0.53788948 -0.06867027 -0.13343620 0.05023384 -0.4124222 0.14546585
## V114 V115 V116 V117 V118 V119
## 16 -0.002893448 0.3607807 -0.03179741 -0.5395298 -0.09937668 0.2623463
## 17 -0.009689331 0.4024983 0.06968880 -0.3447647 -0.06202888 0.5128136
## 18 -0.219610210 0.2596092 0.08492279 -0.2408600 -0.25406265 0.6845303
## 19 -0.272733690 0.2096882 0.02937317 -0.1907062 -0.42234230 0.7405453
## 20 -0.265703200 0.1528912 -0.06367111 -0.1364078 -0.47789001 0.5976582
## 21 -0.139219280 0.1107826 -0.09202957 -0.1172638 -0.40101433 0.4827347
## V120 V121 V122 V123 V124 V125
## 16 0.7113304 0.27961159 0.001560211 0.25109291 0.29890823 -0.03913879
## 17 1.0346889 0.53624916 0.100782390 0.13328171 0.29992104 0.02406311
## 18 1.0982838 0.58554840 0.163515090 -0.01390266 0.16609383 -0.11890602
## 19 1.0828152 0.54117203 0.160310750 0.05171967 -0.00437355 -0.42531204
## 20 1.0239048 0.22397995 0.134063720 0.10375214 -0.28413963 -0.57640457
## 21 0.7621098 -0.01867104 0.166563030 0.15023422 -0.51242065 -0.65453339
## V126 V127 V128 V129 V130 V131 V132
## 16 -0.15289116 -0.5392189 0.3682022 0.5604687 0.5506229 0.46234512 0.36132812
## 17 -0.05968857 -0.6475029 0.4196873 0.6552353 0.9379711 0.61281204 0.18468857
## 18 0.03038979 -0.6403923 0.5849228 0.6591415 1.0959358 0.68452835 0.19828415
## 19 0.12726784 -0.4403114 0.5793724 0.6592941 1.1776581 0.54054451 0.03281403
## 20 0.18429565 -0.0471077 0.3363266 0.5635929 1.0721092 0.29765892 -0.27609444
## 21 0.16077995 0.3607826 0.1079712 0.2827339 0.8989849 0.03273583 -0.18788910
## V133 V134 V135 V136 V137 V138 V139
## 16 -0.4703884 0.3015595 0.2010937 0.34890747 0.21086121 -0.1528931 0.2607803
## 17 -0.4137516 0.5007820 0.3832817 0.34992027 0.12406158 -0.1596889 -0.1475029
## 18 -0.2644539 0.6135139 0.4860954 0.11609459 0.03109169 -0.2696114 -0.6403923
## 19 -0.1588287 0.7603092 0.5517197 -0.00437355 0.02468681 -0.3227329 -0.7403107
## 20 -0.1260185 0.9340630 0.6037521 0.01585960 -0.07640457 -0.2657032 -0.4471092
## 21 -0.2186699 1.0165634 0.6002350 0.23757935 0.04546738 -0.1892204 -0.0392189
## V140 V141 V142 V143 V144 V145
## 16 -0.7817974 0.16046906 -0.49937630 0.03609657 0.7950783 0.8483601
## 17 -0.9803104 -0.04476547 -0.41203117 -0.07968712 0.7746887 0.7725010
## 18 -1.0650768 -0.24085999 -0.35406303 -0.24671745 0.8395328 0.9817963
## 19 -1.0206261 -0.34070778 -0.22234154 -0.53445625 0.7090645 1.1211700
## 20 -0.7636719 -0.43640709 -0.07789040 -0.73484230 0.4301567 1.1452293
## 21 -0.5920296 -0.26726532 0.09898567 -0.64476585 0.1733589 1.0863323
## V146 V147 V148 V149 V150 V151
## 16 1.06906130 0.77859306 0.637657170 0.9833584 0.33085823 -0.02671814
## 17 0.60828018 0.42952919 0.227422710 0.4353123 -0.08218956 -0.30250359
## 18 0.39851379 0.40234566 0.208593370 0.3185940 -0.21085930 -0.16664124
## 19 0.17281151 0.31171989 0.316875460 0.3934364 -0.16648483 0.02218819
## 20 0.01656342 0.10500336 0.229610440 0.3823452 -0.15695190 0.14789200
## 21 -0.05843735 -0.05351639 0.005079269 0.2579651 -0.16296959 0.20453262
## V152 V153 V154 V155 V156 V157
## 16 0.12570190 -0.43577957 -0.8331261 -0.7339058 -0.21742058 -0.7816391
## 17 -0.16781044 -0.68476486 -0.9507809 -0.9159355 -0.37906075 -0.8950024
## 18 -0.09757805 -0.56586075 -0.8340626 -0.8817196 -0.29546738 -0.8582020
## 19 0.09687424 -0.30320549 -0.5198422 -0.6694565 -0.03718758 -0.5488319
## 20 0.26757812 -0.10640717 -0.3866406 -0.4398441 0.32140541 -0.2497692
## 21 0.23421860 -0.08226585 -0.4585152 -0.3460140 0.36960983 -0.2899208
## V158 V159 V160 V161 V162 V163 V164
## 16 -0.30718994 0.1573410 0.05515671 0.3821087 0.6333599 0.91703033 0.61695290
## 17 -0.38046646 0.1532784 0.15867233 0.3565579 0.4078121 0.51374817 0.09093666
## 18 -0.33023643 0.2698460 0.40234566 0.3273430 0.3666382 0.45085907 0.06742287
## 19 -0.20718956 0.2617188 0.58687592 0.2296867 0.2435169 0.09093857 -0.08062744
## 20 -0.00718689 0.2062511 0.67585945 0.3248463 0.2292957 -0.14710999 -0.21367264
## 21 -0.02343750 0.1427345 0.61632729 0.3792152 0.1582813 -0.12171745 -0.06077957
## V165 V166 V167 V168 V169 V170
## 16 0.7117176 1.38562390 1.0110931 -0.05742073 -0.5178909 -0.1209412
## 17 0.2564831 0.51922035 0.5528126 -0.60531235 -1.1012516 -0.3117180
## 18 0.3141422 0.10843658 0.1707783 -1.01671600 -1.3957024 -0.6764889
## 19 0.2755413 -0.09609222 -0.1082058 -1.14843750 -1.2488289 -0.8209400
## 20 0.1223450 -0.07289123 -0.1173420 -1.05734440 -1.0360203 -0.7784386
## 21 0.1727352 0.06648445 0.1189823 -0.86164093 -1.0849228 -0.7984371
## V171 V172 V173 V174 V175 V176
## 16 -0.1339092 -0.4660950 -0.1666393 -0.08039475 0.12578011 0.21570206
## 17 -0.2142200 -0.5988293 -0.2109375 -0.24968910 0.06374741 0.04093742
## 18 -0.4201546 -0.6976547 -0.2576561 -0.20336151 0.23960876 0.02742195
## 19 -0.5032806 -0.5843735 -0.2228146 -0.16773415 0.29718781 0.07187271
## 20 -0.5249977 -0.5328903 -0.1551533 -0.05320358 0.33039093 0.16257668
## 21 -0.4910164 -0.5149231 -0.1007843 -0.12796974 0.20953178 0.11421967
## V177 V178 V179 V180 V181 V182
## 16 -0.22078323 -0.5731258 0.12359619 -0.2949219 0.1958599 -0.1096897
## 17 -0.55101585 -0.8070278 -0.11468697 -0.3703117 -0.2112522 -0.5754681
## 18 -0.51210785 -0.7140636 0.04577828 -0.3154678 -0.5044517 -0.7539864
## 19 -0.30445671 -0.4910946 0.26179314 -0.2834358 -0.6175804 -0.7184391
## 20 -0.08765793 -0.3266411 0.32390594 -0.2535934 -0.4972687 -0.4896889
## 21 -0.11726379 -0.3822651 0.28148270 -0.1966400 -0.1874199 -0.3184357
## V183 V184 V185 V186 V187 V188
## 16 -0.2151585 0.1864052 0.12085915 -0.2991428 -0.2454681 -0.07679748
## 17 -0.5679703 -0.1863270 -0.26968765 -0.8096905 -0.7050018 -0.30031395
## 18 -0.7351551 -0.3064079 -0.31140900 -0.7671108 -0.4891434 -0.26632690
## 19 -0.7907810 -0.2668724 -0.09656334 -0.3864842 -0.1265602 -0.24062729
## 20 -0.5812473 -0.1128921 0.01484489 -0.2807026 -0.1658592 -0.20242310
## 21 -0.4510174 -0.2074223 -0.18203354 -0.5467205 -0.4217186 -0.13328171
## V189 V190 V191 V192 V193 V194 V195
## 16 0.2692184 0.36312485 -0.1226540 0.25008011 0.4158592 0.1753082 -0.4701595
## 17 -0.1335163 -0.11702919 -0.5834370 0.04843903 0.3225002 0.2020321 -0.2867184
## 18 -0.2958603 -0.41156387 -0.7179718 -0.14296722 0.2180481 0.3510132 -0.1589031
## 19 -0.4294567 -0.42484093 -0.5457058 -0.15093803 0.2024193 0.4028091 -0.2270317
## 20 -0.3414059 -0.17914009 -0.2060928 -0.10359383 0.2714806 0.3728123 -0.3724976
## 21 -0.1660156 -0.06351471 -0.1697655 -0.11788940 0.3925800 0.4190617 -0.2622662
## V196 V197 V198 V199 V200 V201
## 16 0.07890701 0.6096096 1.0221043 0.47328186 0.4857006100 0.58547020
## 17 0.26742172 0.7765636 0.7503109 0.12499809 -0.0790615080 0.18398285
## 18 0.19734573 0.6723423 0.4016399 0.03335762 -0.1675758400 0.01664162
## 19 -0.10312271 0.3846874 0.1810169 0.13718987 0.0006217956 -0.07320786
## 20 -0.36289024 0.2610950 0.2867947 0.21664047 0.0688266750 -0.10890579
## 21 -0.30367279 0.4079666 0.4595299 0.26703262 0.0092182159 -0.18976593
## V202 V203 V204 V205 V206 V207
## 16 0.9256229 0.37484550 0.67757797 0.61710930 0.73031235 0.3635922
## 17 0.7967205 -0.08593941 0.05968666 0.24624825 0.32953262 -0.1304684
## 18 0.7534370 -0.25172234 -0.23921776 0.06679535 -0.01148605 -0.4939041
## 19 0.6614075 -0.29195595 -0.38718796 -0.09132957 -0.30593872 -0.6345291
## 20 0.5396099 -0.11484337 -0.43234444 -0.18851852 -0.26968765 -0.6574974
## 21 0.4102344 0.05648231 -0.55289078 -0.37616920 -0.16718864 -0.5597649
## V208 V209 V210 V211 V212 V213
## 16 0.3026562 0.7621098 0.22335625 0.12828064 0.67695236 0.49046707
## 17 -0.4175777 0.2878094 0.03155899 0.04874802 0.57343674 0.14898491
## 18 -0.8276558 -0.0701580 -0.08461189 0.03085899 0.35367203 -0.01460838
## 19 -0.8168735 -0.1740627 -0.01273155 0.22218895 0.28437424 -0.01195717
## 20 -0.8103924 -0.3301563 0.02929497 0.25163841 0.17507553 -0.03765678
## 21 -0.6874218 -0.2782841 -0.07547188 0.06828117 0.09671974 -0.04101562
## V214 V215 V216 V217 V218 V219
## 16 0.2331219 0.1260948 0.5738277 0.1471119 -0.01969147 -0.3151588
## 17 -0.1932812 -0.5146885 -0.0465622 -0.3875027 -0.53296852 -0.7954674
## 18 -0.3840637 -0.9417210 -0.4767170 -0.7269554 -0.49648476 -0.8614063
## 19 -0.1685925 -0.6757069 -0.2846870 -0.5325794 -0.13094139 -0.6257820
## 20 -0.1003914 -0.3348427 -0.1773453 -0.4710197 -0.06593704 -0.5387478
## 21 -0.1497650 -0.3322658 -0.4166412 -0.5811691 -0.15468788 -0.5660152
## V220 V221 V222 V223 V224 V225 V226
## 16 0.11265945 0.7871094 0.18210983 0.6070309 0.3407040 0.5929680 1.3456230
## 17 -0.18507957 0.4715614 0.07031059 0.3962479 0.1459370 0.2514858 1.1279716
## 18 -0.29140663 0.3098431 0.12163925 0.4771080 0.1499233 0.1228924 0.9871864
## 19 -0.14812088 0.4196873 0.26851654 0.6059380 0.3931236 0.2492924 1.0926571
## 20 -0.05539131 0.5448437 0.30429459 0.7591400 0.8350773 0.6285915 1.2883587
## 21 -0.22867203 0.4192162 0.22702980 0.7982826 1.0704689 0.8977356 1.4464855
## V227 V228 V229 V230 V231 V232 V233
## 16 0.3835945 0.4400787 0.23710823 -0.4234409 -0.6201553 0.49265671 0.7646084
## 17 0.6715622 0.4509373 -0.26375008 -0.9267159 -0.9779682 0.07742119 0.4590607
## 18 0.9320297 0.3357830 -0.43195343 -0.8952370 -0.8676548 -0.09015846 0.2360935
## 19 1.1880436 0.3053131 -0.28133011 -0.6159382 -0.6045322 0.02437592 0.3471889
## 20 1.4201565 0.3601570 -0.25977135 -0.5934372 -0.7362499 -0.07538986 0.3373451
## 21 1.5002327 0.6358604 -0.03116989 -0.4821873 -0.7747669 -0.07367325 0.3992157
## V234 V235 V236 V237 V238 V239
## 16 0.36210823 -0.10046959 -0.22554779 -0.773281100 -0.6468754 -0.34015465
## 17 -0.03343773 -0.38250160 -0.48406219 -1.066015200 -1.0020313 -0.64093590
## 18 -0.10336113 -0.26414108 -0.39007568 -0.840858460 -0.8053150 -0.43421936
## 19 0.11226654 0.04843903 -0.01687622 -0.300706860 -0.3085938 -0.03320503
## 20 0.11304474 0.13039017 0.27382660 0.001092911 -0.1491413 0.03390503
## 21 0.13828087 0.21203232 0.44046974 0.072734833 -0.1610146 -0.08726502
## V240 V241 V242 V243 V244 V245
## 16 -0.4949207 -0.3678913 0.09031105 -0.381408690 0.08140755 0.21586037
## 17 -0.9828110 -1.0800018 -0.52671814 -0.732967380 -0.08757973 0.14530945
## 18 -1.0229645 -1.3569527 -0.57023621 -0.492654800 -0.01640701 0.04109383
## 19 -0.8709374 -1.1675797 -0.16843987 -0.007030487 0.13937759 0.02593803
## 20 -0.9298458 -1.0422707 0.01156044 0.127500530 0.16585922 0.16359520
## 21 -1.0541401 -0.9186726 0.12406158 0.223983760 0.25757790 0.41796684
## V246 V247 V248 V249 V250 V251 V252
## 16 0.34710693 0.64578056 -0.1067963 -0.1845303 0.2131214 -0.3064060 0.5575790
## 17 0.01656151 0.21499825 -0.4565620 -0.6985149 -0.4407787 -1.0421886 0.1659374
## 18 -0.32461166 0.03210831 -0.5163288 -0.7821102 -0.6190624 -1.1242218 0.2407837
## 19 -0.33648300 0.09718704 -0.3043766 -0.5294571 -0.5198422 -0.8394566 0.4990635
## 20 -0.15070343 0.13164139 -0.3011723 -0.3976574 -0.4916420 -0.5148430 0.6739044
## 21 -0.13921928 0.08828354 -0.2282810 -0.3597641 -0.5772648 -0.5585155 0.5121098
## V253 V254 V255 V256 V257 V258
## 16 0.6883602 0.52280807 0.1173420 -0.3673420 0.04461098 0.555856700
## 17 0.3287487 0.31078148 -0.3504696 -0.6413288 0.05531120 0.230312350
## 18 0.2892971 0.30726242 -0.5126534 -0.7239075 -0.05515671 0.002887726
## 19 0.3824196 0.26655960 -0.6420307 -0.6206226 -0.09406090 -0.083982468
## 20 0.3589802 0.09031296 -0.9037476 -0.6853905 -0.25265503 -0.118204120
## 21 0.2575779 -0.02968788 -1.1035175 -0.8199196 -0.34328461 -0.152969360
## V259 V260 V261 V262 V263 V264
## 16 0.3945312 0.385702130 0.74921799 0.65187263 -0.1651535 -0.5974217
## 17 0.1187477 -0.029062271 0.26773453 0.37172127 -0.1821861 -0.7328110
## 18 0.1196079 -0.096326828 0.13914108 0.29468536 -0.3179703 -0.9342156
## 19 0.1959381 0.005624771 0.04554367 0.17140770 -0.4232063 -0.9621849
## 20 0.2053909 -0.008672714 -0.13515854 -0.00164032 -0.3598423 -0.9435921
## 21 0.1345310 -0.067031860 -0.07351494 0.04273605 -0.3147659 -1.0328903
## V265 V266 V267 V268 V269 V270
## 16 -0.2053909 0.006557465 -0.4001579 -0.07859039 0.63336182 0.638355260
## 17 -0.3000011 -0.175468440 -0.7092209 -0.43757820 0.20406151 -0.005937576
## 18 -0.4544544 -0.065237045 -0.6789055 -0.42765617 -0.06265831 -0.500860210
## 19 -0.4400787 0.131561280 -0.5432816 -0.34187317 -0.23531342 -0.660232540
## 20 -0.3622704 0.125312810 -0.7737465 -0.52038956 -0.39390564 -0.603204730
## 21 -0.2936726 0.071561813 -1.0147648 -0.67117310 -0.42453384 -0.462970730
## V271 V272 V273 V274 V275 V276 V277
## 16 0.84953117 0.7444515 1.0642185 0.7081222 0.2785950 -0.02117157 0.24335861
## 17 0.60374832 0.5446892 0.6114845 0.1242218 -0.4334373 -0.40781212 0.29999924
## 18 0.22085762 0.2886734 0.1203918 -0.2240620 -0.7629700 -0.70046616 0.47054482
## 19 0.04843903 0.3018723 -0.2057056 -0.5110931 -0.8394566 -0.84718704 0.55867004
## 20 0.03413963 0.3513260 -0.3414059 -0.8741417 -0.7823429 -0.87984467 0.35023117
## 21 0.04578209 0.2629700 -0.3422661 -0.8860149 -0.5047665 -0.80289078 0.06883049
## V278 V279 V280 V281 V282 V283
## 16 0.1703110 -0.76391029 -0.133592610 0.282110210 0.59835815 1.1932812
## 17 0.5270328 -0.33671761 -0.086328506 0.086561203 0.25656319 0.8962479
## 18 0.8060112 0.06734657 0.007341385 0.069843292 0.26788902 1.0121078
## 19 0.7203102 0.02672005 -0.040622711 0.004686356 0.26101685 1.0909367
## 20 0.3765640 -0.29874802 -0.186639790 -0.167654040 0.15929413 1.1091404
## 21 0.1015606 -0.49976540 -0.241174700 -0.243284230 0.02453041 1.0557804
## V284 V285 V286 V287 V288 V289 V290
## 16 1.634453 1.1504688 0.60812187 0.07484436 -0.01992035 0.5021095 -0.02719116
## 17 1.203436 0.5914841 0.08797073 -0.57593536 -0.61906242 0.2437477 -0.38671684
## 18 1.213673 0.6216392 0.12218857 -0.71796989 -0.85421753 0.2405472 -0.51648712
## 19 1.271874 0.6980419 0.38390732 -0.40945625 -0.71343613 0.5386696 -0.34093857
## 20 1.145077 0.5698433 0.44335938 -0.06484222 -0.51359558 0.7352295 -0.25093651
## 21 1.000469 0.4539852 0.40023613 0.03898239 -0.39163971 0.7938290 -0.18718910
## V291 V292 V293 V294 V295 V296
## 16 -0.4151573 -0.4085922 0.10335922 0.79585648 0.913280490 0.58570290
## 17 -0.6767178 -0.7700787 -0.34844017 0.48280907 0.487499240 0.18218803
## 18 -0.6276531 -0.9126549 -0.44265747 0.25413895 0.223356250 0.06117058
## 19 -0.4332790 -0.8231239 -0.18906403 0.14476585 0.028438568 0.13812065
## 20 -0.3974953 -0.6878910 -0.05265427 0.17429543 0.006639481 0.19507599
## 21 -0.3435135 -0.6411705 -0.18328476 0.06202888 -0.096719742 0.13296890
## V297 V298 V299 V300 V301 V302
## 16 0.23546791 -0.09062576 0.38109589 0.18008041 0.21461105 -0.06093979
## 17 -0.09226608 0.02172089 0.21531296 -0.11156464 0.08499908 -0.05421829
## 18 -0.11460876 0.03718567 0.08827972 -0.13046646 -0.08820343 0.00601387
## 19 -0.03070831 -0.10734367 -0.05195618 -0.08093643 -0.04382896 0.07156181
## 20 -0.05515671 -0.17663956 0.10515785 -0.05109405 0.30648041 0.20906258
## 21 -0.10726547 -0.19726753 0.10523415 -0.25164032 0.42882919 0.30406380
## V303 V304 V305 V306 V307 V308 V309
## 16 0.29359055 0.29015541 0.7721100 0.75460815 0.70828056 0.8419533 0.4429684
## 17 0.05327988 0.02116966 0.5815601 0.66781044 0.70374870 0.7959385 0.7502346
## 18 -0.33140373 -0.27640724 0.3110924 0.45788956 0.54585838 0.7011719 1.0791397
## 19 -0.54827881 -0.31187248 0.1621876 0.23226547 0.25843811 0.5168724 0.8530445
## 20 -0.41749954 -0.19038963 0.2785950 0.20304489 0.10039139 0.5563278 0.6723423
## 21 -0.29351616 -0.10367203 0.2454662 0.03703117 0.09703064 0.3904686 0.3277340
## V310 V311 V312 V313 V314 V315 V316
## 16 -0.1281261 0.1448460 0.5113296 1.208361 0.3628101 -0.19390678 0.65640640
## 17 0.2329693 0.4465618 0.7159386 1.297499 0.4120312 -0.48546982 0.19867325
## 18 0.6484375 1.0332794 1.1445332 1.378048 0.7472629 -0.40265274 0.09484482
## 19 0.4576569 1.0342960 1.1665630 1.243673 0.8503094 -0.17453003 0.38437843
## 20 0.3858585 1.1189060 1.2851543 1.128979 1.1790619 0.08875275 0.61460876
## 21 0.2927342 1.0514832 1.3158588 1.017578 1.2753105 -0.02476311 0.25757790
## V317 V318 V319 V320 V321 V322 V323
## 16 0.59710884 1.1371059 0.6007824 0.7444515 -0.1945305 0.01437187 0.13234520
## 17 0.35156059 1.0865612 0.7849960 0.7121887 0.1389847 0.21921730 0.05531120
## 18 0.06984329 0.8091392 0.8708592 0.7199230 0.5978908 0.38843536 -0.19922256
## 19 0.09718895 0.4385166 0.5584393 0.5481224 0.6417923 0.31515694 -0.35320473
## 20 0.28609276 0.3030453 0.3666420 0.4713268 0.6698418 0.49085999 -0.07734299
## 21 0.11046600 0.1432781 0.2645302 0.3579712 0.5952339 0.47398567 -0.16976547
## V324 V325 V326 V327 V328
## 16 0.43882942 -0.3291397 0.7328091 0.8723450 0.4501572
## 17 0.23093796 -0.7037506 0.2107830 0.5207806 0.2324200
## 18 -0.06046677 -1.2319546 -0.6264858 -0.1476555 -0.1301556
## 19 -0.04093361 -1.1313286 -0.6496925 -0.2182808 -0.0706253
## 20 0.04640579 -0.7122688 -0.5821876 -0.3124962 -0.1203918
## 21 -0.22789192 -0.8699207 -0.7609367 -0.5010166 -0.3286724
t(SSTdata) %>% dim()
## [1] 328 2261
# construct the EOFs
Z <- t(SSTdata) # data matrix, rows are times, columns are spatial sampling locations.
spat_mean <- apply(SSTdata, 1, mean) # average across all time.
nT <- ncol(SSTdata) # no of time points
# outer command makes the spatial mean into a matrix
Zspat_detrend <- Z - outer(rep(1, nT), spat_mean) # detrend data (center to spatial mean)
Zt <- 1/sqrt(nT-1)*Zspat_detrend # normalize
E <- svd(Zt) # SVD
# number of EOFs to estimate in the data
n <- 10
# Projected onto the
TS <- Zspat_detrend %*% E$v[, 1:n]
summary(colMeans(TS))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.091e-16 -8.031e-17 3.201e-16 3.113e-16 6.916e-16 8.885e-16
This verifies that the time series have a zero mean.
tau <- 6
nT <- nrow(TS)
TStplustau <- TS[-(1:tau), ] # TS with first tau time pts removed
TSt <- TS[-((nT-5):nT), ] # TS with last tau time pts removed
Cov0 <- crossprod(TS)/nT # cross product
Covtau <- crossprod(TStplustau, TSt) / (nT - tau)
C0inv <- solve(Cov0)
Mest <- Covtau %*% C0inv # estimated M, the dynamics
Ceta <- Cov0 - Covtau %*% C0inv %*% t(Covtau) # Covariance matrix
image(Mest)
image(Ceta)
The forecast is given by \(\hat\mu +
\PhiM^2\alpha_t\).
SSTlonlat$pred <- NA
alpha_forecast <- Mest %*% TS[328, ] # propigate alpha first,
idx <- which(SSTlonlat$mask == 0)
SSTlonlat$curr[idx] <- as.numeric(E$v[, 1:n] %*% TS[328, ] + spat_mean)
SSTlonlat$pred[idx] <- as.numeric(E$v[, 1:n] %*% alpha_forecast + spat_mean) # add to spatial mean for original scale.
# for plotting
SSTlonlat$obs1[idx] <- SSTdata[, 328]
SSTlonlat$obs2[idx] <- SST_Oct97
C <- Mest %*% Cov0 %*% t(Mest) + Ceta # 6 month ahead variances
SSTlonlat$predse[idx] <- sqrt(diag(E$v[, 1:n] %*% C %*% t(E$v[, 1:n])))
SSTlonlat %>% ggplot(aes(lon, lat, fill = pred)) +
geom_tile()
See the appendix for proper visualization.
DSTM_EM is provided with the package STRbook, and runs
an EM algorithm that does ML estimation in state-space models.
DSTM_Results <- DSTM_EM(Z = SSTdata,
Cov0 = Cov0,
muinit = matrix(0, n, 1),
M = Mest,
Ceta = Ceta,
sigma2_eps = 0.1,
H = H <- E$v[, 1:n],
itermax = 10,
tol = 1)
The object DSTM_Results has the setimated parameters,
and complete-data negative log-likelihood. Estimates of \(\alpha_{1,t}\), \(\alpha_{2,t}\) and \(\alpha_{3,t}\)
par(mfrow = c(1,3))
for(i in 1:3) {
plot(DSTM_Results$alpha_smooth[i, ], type = "l",
xlab = "t", ylab = bquote(alpha[.(i)]))
lines(TS[, i], lty = "dashed", col = 'red')
}
image(as(DSTM_Results$Mest, "Matrix"))
image(as(DSTM_Results$Mest %^% 6, "Matrix"))
image(as(Mest, "Matrix"))
alpha <- DSTM_Results$alpha_smooth[, nT]
P <- DSTM_Results$Cov0
for(t in 1:6) {
alpha <- DSTM_Results$Mest %*% alpha
P <- DSTM_Results$Mest %*% P %*% t(DSTM_Results$Mest) + DSTM_Results$Ceta
}